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Abstract 

We use numerical simulations to address locomotion at zero Reynolds number in viscoelastic 
(Giesekus) fluids. The swimmers are assumed to be spherical, to self-propel using tangential 
surface deformation, and the computations are implemented using a finite element method. The 
emphasis of the study is on the change of the swimming kinematics, energetics, and flow disturbance 
from Newtonian to viscoelastic, and on the distinction between pusher and puller swimmers. In 
all cases, the viscoelastic swimming speed is below the Newtonian one, with a minimum obtained 
for intermediate values of the Weissenberg number. We. An analysis of the flow field places the 
origin of this swimming degradation in non-Newtonian elongational stresses. The power required 
for swimming is also systematically below the Newtonian power, and always a decreasing function 
of We. A detail energetic balance of the swimming problem points at the polymeric part of the 
stress as the primary VFe-decreasing energetic contribution, while the contributions of the work 
done by the swimmer from the solvent remain essentially VFe-independent. In addition, we observe 
negative values of the polymeric power density in some flow regions, indicating positive elastic work 
by the polymers on the fluid. The hydrodynamic efficiency, defined as the ratio of the useful to 
total rate of work, is always above the Newtonian case, with a maximum relative value obtained at 
intermediate Weissenberg numbers. Finally, the presence of polymeric stresses leads to an increase 
of the rate of decay of the flow velocity in the fluid, and a decrease of the magnitude of the stresslet 
governing the magnitude of the effective bulk stress in the fluid. 
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I. INTRODUCTION 



The physical mechanisms of microorganism locomotion accompany a variety of natural 
phenomena, including human spermatozoa approaching the ovum in the mammalian female 
reproductive tract [Ij, algal blooms moving to nutrient rich environment in maritime regions 
[2l|3j, biofilm formation [4J, and paramecia cells escaping from their predators [5j. Much work 
has been done to shed light on the different physical mechanisms at play in the locomotion 
of microorganisms, including classical kinematics studies P-Ej, nutrients uptake ITT] , 
collective behavior [T2HT6] and hydrodynamic interactions [T7H19]. 

Most past work has however been limited to the study of locomotion in Newtonian flu- 
ids, whereas in many biological instances microorganisms swim in complex polymeric fluids. 
Several groups have recently started to address the effect of viscoelasticity on the loco- 
motory features and hydrodynamic performance of swimming microorganisms. Lauga ^Oj 
investigated the kinematics and energetics of Taylor's waving sheet in various model non- 
linear polymeric fluids and found that, for a given waving stroke, viscoelasticity impeded 
the locomotion of the swimming sheet. Similarly, Fu et al. [21j derived analytically that 
the swimming speed of an inflnite fllament with a prescribed waveform will decrease in 
Oldroyd-B fluid compared to the Newtonian case. Teran et al. |22j numerically studied 
the influence of viscoelasticity on the dynamics of a flnite-length Taylor's waving sheet and 
showed that the swimming speed and efficiency can actually increase within a certain range 
of polymer relaxation time and sheet waveforms. Recently, Shen et al. [23j experimentally 
demonstrated that fluid elasticity hindered the self-propulsion of the nematode Caenorhab- 
ditis elegans^ primarily due to polymeric stretching near hyperbolic points in the viscoelastic 
flow. 

In our previous work [24j, we carried out three-dimensional axisymmetric simulations 
to study the effect of viscoelasticity on the locomotion speed and energy consumption of 
a spherical squirmer swimming by tangential surface deformations, as a model of ciliated 
protozoa. A non-linear viscoelastic model (Giesekus) was used to describe the polymer dy- 
namics. We showed the dependency of swimming speed, power consumption, and swimming 
efficiency on the fluid elasticity. This work was however restricted to the most energetically 
efficient swimming gait, one associated with no vortices generated around the cell in the 
Newtonian case (a potential flow squirmer in the Newtonian limit). 
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In the current paper we numerically investigate the hydrodynamic performance of a 
squirmer utilizing other swimming gaits in a viscoelastic fluid. The gaits correspond to 
pusher and puller swimmers, a natural distinction for biological organisms. We assume 
the gait to be steady and thus ignore, as a flrst model, the oscillatory flow fleld around a 
swimming cell. As recently shown by Guasto et al. and Drescher et al. [26], the flow fleld 
generated by some microorganisms such as Chlamydomonas reinhardtii cannot be accurately 
described by a steady analysis, and in fact, the flow generated by many swimmers is usually 
time-periodic, with large fluctuations around the mean. However, the simple steady model 
we utilize here has been used to understand several fundamental processes related to the 
physics of swimming microorganisms, such as nutrient uptake [27j, locomotion in stratifled 
fluid [28j , biomixing [29j and the collective behavior of microorganisms p!4l l30 j . Note that 
in the case of viscoelastic fluids like those studied here, the presence of an additional time 
scale (stroke frequency) leads to the deflnition of a Deborah number, different from the 
Weissenberg number deflned below based on shear rates, and possibly leading to a rich 
dynamical behavior. 

The paper is organized as follows. We flrst investigate in detail how fluid elasticity affects 
the swimming speed of the spherical model cells applying different locomotive gaits and the 
resulting polymeric dynamics in the flow fleld. We then analytically derive a decomposition 
of the power associated to cell locomotion by surface deformation in the viscoelastic fluid, 
extending the analysis presented for the Newtonian fluid in Ref . [7] . We further compute the 
change in the swimming speed of the squirmer in the viscoelastic fluid under const ant- power 
conditions, as well as the hydrodynamic efficiency. Finally, in order to quantify the influence 
of the squirming motion on nearby swimmers as well as on the bulk stress of the surrounding 
fluid, we analyze the velocity decay around the model swimming cells and calculate their 
effective stresslets. 

II. MATHEMATICAL MODEL 
A. Squirmer Model 

In this so-called envelope approach, flrst proposed by Blake [3T], one assumes the squirmer 
imposes a non-zero tangential velocity on its surface, U5, in the co-moving frame. Here we 
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adopt the concise formulation introduced in Ref. Jlj 



/ N 2 _ / /e-r\ /e-rr \ 
'^^(''^ = J2^^^^-^nl ( e , 1 

n>l ^ ^ 

where e is the orientation vector of the squirmer, is the nth mode of the surface squirming 
velocity [31j, Pn is the nth Legendre polynomial, r is the position vector, and r = |r|. In 
a Newtonian fluid, the swimming speed of the squirmer is U^ew = 25i/3 [31j and thus 
only dictated by the flrst mode. As in many previous studies [TOl [HI [TTl [T^j, we assume 
5^ = for n > 3. The tangential velocity on the sphere in the co-moving frame is therefore 
expressed as ue{9) = Bis'mO + (S2/2)sin20, where 9 = arccos(e • r/r). We introduce an 
additional parameter a, representing the ratio of the second to the flrst squirming mode, 
thus a = B2/B1. When a is positive, the swimmer get impetus from its front part, and 
is termed a puller (swimmers in this category include the alga genus Chlamydomonas) . As 
a difference, when a is negative, the thrust comes from the rear part of the body and the 
swimmer is a pusher (swimmers in this category include all peritrichous bacteria such as 
Escherichia coli). If not otherwise stated, here we will present results obtained with a = —5 
(simply denoted as pusher) , a = 5 (puller) and a = (neutral squirmer inducing a potential 
velocity fleld in the Newtonian case). 



B. Polymeric Dynamics 

For steady incompressible low- Reynolds number flow in a viscoelastic fluid, the momen- 
tum and continuity equation are written as 

- Vp + V-r = 0, (2) 
V • u = 0. (3) 

The velocity is scaled by Bi ^ related by the expressions above to the magnitude of the 
velocity of the boundary and to the swimming speed while lengths are scaled by the 
diameter of the spherical swimmer time is scaled by D/Bi^ and pressure and stresses are 
scaled by jiBi/D^ where ji is the coeflicient of viscosity. 

Following classical modeling approaches f32] - [34j . the deviatoric stress r is split into two 
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components, the viscous solvent stress (r^) and the polymeric stress (r^), so r = + r^. 
The stress is governed by Newton's law of viscosity, thus is formulated as 

t' = (3 (Vu + Vu^) , (4) 

where /5 < 1 represents the ratio of the solvent viscosity, /z^, to the total zero shear rate vis- 
cosity, /X, of the solution. To close the model, a transport equation for the polymeric stress 
is required. Here we adopt the nonlinear Giesekus model [l35J, which, in addition to 
displaying shear-thinning material properties (viscosity and normal stress differences), pro- 
vides two important features, namely saturation of polymer elongation, and a non-negative 
entropy production during the time evolution of the polymers (see details in Ref. [361438] ). 
Violation of these two properties may cause non-physical flow behavior as well as numerical 
difliculties. The nondimensional constitutive equation for this model can be written as 

rP^WeTP + -^—-f{TP-T^)= (l-/5)(Vu + Vu^), (5) 

V 

where A denotes the upper-convected derivative, deflned for a tensor A by 

V <9A 

A = — + u • VA - Vu^ . A - A . Vu. (6) 

dt ^ ^ 

In the constitutive equation for the polymeric stress, Eq. ([5]), We denotes the Weissenberg 
number, deflned as We = XBi/D^ where A is the polymer relaxation time. The so-called 
mobility factor is introduced in the nonlinear stress term to represent an anisotropic 
hydrodynamic drag on the polymer molecules [32] and limits the extensional viscosity of the 
fluid. From thermodynamics considerations, the mobility factor a^n must be in the to 1/2 
range [32[ |39], and in this paper we assume = 0.2. 



III. NUMERICAL METHOD 

A numerically stable and accurate flnite element model is built, based on the formulation 
denoted Discrete Elastic- Viscous Split Stress (DEVSS-G) [40l[41^. The governing equations 
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are rewritten as 



V . fiaiVu + Vu^) - + V . - V . (/X, - /3) (G + G^) = 0, (7) 

where an additional tensor G = Vu is introduced as an independent interpolation of the 
velocity gradient tensor Vu and an additional elliptic term V-/Xa(Vu+Vu^) — V-/ia(G+G^) 
is added into the momentum equation [42j. In our computations, /x^ is chosen to be 1 as in 
Liu et al. [41j, and the velocity gradient term in equation for the polymeric stress, Eq. ([5]), 
is approximated by G. 

A Galerkin method is used to discretize the momentum equations, continuity equation, 
and the equation for the additional unknown G. We use quadratic elements for u and 
linear elements for both p and G. To improve the numerical stability, the streamline- 
upwind/Petrov-Galerkin(SUPG)[43j method is used to discretize the constitutive equation, 
Eq. ([5]). Finally, the weak form for the constitutive equations is written as 

|S + -^u . VS, TP+We{u . Vr^ - G^ . - . G + -^(r^ • r^)) 
^ Uc 1 — P 

- (l-/5)(Vu + Vu^)}=0, (8) 

where S denotes the test function for r^, /i is the characteristic length-scale of the element 
and t/c is the magnitude of the local characteristic velocity. In our case, we choose the norm 
of u as the value of Uc- The finite-element framework for our implementation is provided 
by the commercial software COMSOL. 

We perform three-dimensional axisymmetric simulations on a two-dimensional mesh. Tri- 
angle elements are used for all the simulations and the number of elements are dependent 
on the swimming gait. We generate sufficiently refined mesh in regions with high magnitude 
of polymeric stress to overcome numerical instabilities |^ HS] , In all cases, high resolution 
near the cell body is required in order to resolve the thin stress boundary layers. Fine meshes 
are also necessary to capture the high polymer elongation rates. Typical mesh size at the 
cell boundary is 5E10~^ of the reference (swimmer) length. The total number of elements 
used for the results reported here ranges from 150, 000 to 250, 000. 

The numerical implementation is first validated against the theoretical results for the 
swimming speed and power of squirmers of different gaits in Newtonian ffuids [31j. For 
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FIG. 1. (Color online) Vector plot of flow fields generated by pusher and puller squirmers in a 
Newtonian fluid. Left: co-moving frame. Right: laboratory frame. On each panel, the data for 
the pusher are displayed on the left, and those for the puller on the right. The large arrow (blue 
online) indicates the swimming direction and the background color scheme indicates the velocity 
magnitudes. 

viscoelastic fluids, we validate the present numerical approach against the simulations in 
Ref. [46] of a sphere sedimenting in a tube filled with an Oldroyd-B fluid (the mobility 
factor a^n equals to zero) and against the results in Ref. [47j where the authors numerically 
study the negative wake of a sphere sedimenting in a tube filled with a Giesekus fluid. 



IV. SVS^IMMING SPEED IN VISCOELASTIC FLUIDS 

To illustrate the difference between the swimming modes considered, the flow fields gen- 
erated by a pusher and a puller in a Newtonian fluid are shown in Fig. [l} where data for 
the pusher are displayed on the left of each panel, and those for the puller on the right. 
Vectors display the in-plane velocity magnitude and direction while the background color 
indicates velocity magnitudes. We see an axisymmetric vortex ring generated in the front of 
the pusher and behind the puller; as a difference the neutral swimmer generates a symmetric 
velocity field with no vorticity |24j. 

To address the dynamics in a viscoelastic fluid, we first examine the dependency of the 
squirmer swimming speed, [/, on the Weissenberg number. We. The swimming speed is 
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FIG. 2. (Color online) Nondimensional swimming speed, ?7, of the swimmer as a function of the 
Weissenberg number, VFe, for pusher a — (squares, green online), puller a — ^ (triangles, blue 
online) and neutral squirmer a — {) (circles, red online). The swimming speed is scaled by its 
corresponding value in the Newtonian fluid, ?7New 

computed as in Ref. [24j. Simulations are performed in the co-moving frame with diflFerent 
values of the inflow velocity, and the total force on the body is computed. The swimming 
speed is then estimated by interpolating to the inflow velocity yielding zero net force on the 
swimmer. A new simulation is then performed to ensure that the value of the predicted 
swimming speed used as free-stream velocity indeed gives values of the total force below 
a given tolerance. The numerical results are shown in Fig. [2| for three diflFerent swimming 
gaits, where the swimming speed is scaled with the speed of the Newtonian squirmer, t/New = 
2Si/3. 

Considering the data in Fig. [2} we flrst observe that viscoelasticity clearly distinguishes 
between the three swimming gaits in terms of [/, while they produce the same swimming 
speed in the Newtonian fluid. The neutral gait, a = 0, generally gives the highest value 
of U for most We numbers while pushers generally have the lowest velocity. Second, for 
all the gaits and We numbers considered here, the swimming speed of swimmers in the 
viscoelastic fluid is lower than that of swimmers in the Newtonian fluid. Finally, for all three 
gaits, U systematically displays a minimum in the range of We considered. The swimming 
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FIG. 3. (Color online) Distribution of the t^z component of the polymeric stress for the neutral 
squirmer and three values of the Weissenberg number, We — 0.1, 2, and 6. (Left: We — 0.1 and 2. 
Right: H^e = 6 and 2.) 

speed initially decreases as We increases from 0, while it eventually always increase for the 
largest values of We. The value of the Weissenberg number corresponding to the minimum 
swimming speed varies with the gait: it increases from 0.5 to 4 when moving from puller to 
neutral to pusher. 

Different swimming gaits generate different polymer dynamics around the swimmer and 
this, in turn, influences the swimming speed as shown above. To gain insight into the 
numerical results, we start by showing in Fig. |3]the distribution of axial polymeric stress, 
Tzz^ for the neutral swimmer [a = 0) in the cases We = 0.1 and We = 2. The value 
We = 0.1 is chosen since it is close to the Newtonian case {U = 0.987), while U is close to 
the minimum swimming speed when We = 2. As seen in the figure, there exists a significant 
difference in the Tzz distribution behind the swimmer. In the case We = 2, the magnitude 
of Tzz is much higher, indicating a higher extent of polymer stretching in this region. As 
suggested in Fig.[l} an elongational flow [32j is generated aft the cell body. In the viscoelastic 
fluid such flow would be responsible for polymer streching which in return increases the local 
elongational viscosity [48l[49j, yielding a strong elastic resistance against locomotion. We 
further see in Fig.[3]that, upon increasing the polymer relaxation time to We = 6, the region 
of largest elongation becomes narrower and closer to the pole, leading to a reduced elastic 
resistance and consistent with the results of Fig. [2j 

The polymer contribution for the pusher swimmer is displayed in Fig. [i] (left). In this 
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case, the kinematics of surface deformation draws fluid from the side and pushes it towards 
the two poles, resulting in strong elongational flows at the poles (see Fig.[T]). Such flow aligns 
the polymer chains near the poles in the swimming direction, leading to high elongational 
viscosities in the front of and behind the squirmer. Comparing Fig. [3] and Fig. [4j we see 
that the magnitude of r^^ on the back is much larger in the case of pushers. This is further 
further illustrated in Fig. [5] where we display the value of the normal polymeric stress along 
the symmetry axis. It is noteworthy that for the pusher swimmer, signiflcant r^z is evident 
also in the region ahead of the body. In this case, the stretched polymers might contribute 
an elastic driving force on the swimmer. However, the magnitude of r^z at the front is lower 
than that behind, and thus the net effect is that of an additional drag. Note that for the 
neutral squirmer, the region of polymer elongation becomes narrower at large We and we 
observe a recovery of the swimming speed. 

We then consider the puller swimmer, which is about 20% faster than the pusher (as seen 
in Fig. [2]). The polymeric stresses in this case are shown in Figs, [i] (right), [5] and [6| The 
results in Fig. |6]are displayed for We = 0.5, corresponding to the minimal swimming speed 
for the viscoelastic puller (and thus the maximum influence of viscoelastic stresses). In the 
case of pullers, the inward velocity imposed by the surface deformation on the front and 
back of the body takes fluid away from the back of the swimmer, and thus stretched fluid is 




FIG. 4. (Color online) Distribution of the Tzz component of the polymeric stress for the pusher 
(left) and puller (right) swimming in the viscoelastic fluid. For both gaits, We = 0.1 and We = 2 
are chosen for comparison. The inset plots show regions with high value of Tzz- 
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FIG. 5. (Color online) Distribution of the r^z component of the polymeric stress along the symmetry 
axis. Shaded circles and arrows indicate the squirmer and its swimming direction. Data for the 
pusher are shown in the top whereas the neutral squirmer and the puller are reported at the bottom. 
The solid and short-dashed lines correspond to We = 0.1, the long-dashed and dot-dashed lines to 
We ^2. 



removed from the region behind the body. As shown in Fig. [5} r^^ is smaller for the puller 
than for the neutral squirmer in the region immediately behind the body. One may thus 
expect the puller may swim faster than the neutral squirmer. However, large magnitude of 
and Trz are observed in conjunction to the vortex ring, with maximum values attained 
near the swimmer surface. Both components contribute to a polymeric (elastic) force applied 
on the swimmer as indicated by the grey arrows in Fig. [6) We observe an elongational flow 
generated on top of the vortex ring, and polymer chains stretched by such flow result in a 
strong force opposing the swimmer motion, possibly explaining why the puller swimmer is 
not faster than the neutral squirmer despite displaying a signiflcant reduction in polymer 
stretching behind its body. 

To summarize, we observe elongational flows generated by all three swimmers around 
their bodies (neutral squirmer, pusher, and puller). The strength, orientation and position 
of the elongational flow is dependent on the swimming gaits. Such flows yield increasing 
elongational viscosities which serve as an additional elastic force, possibly positive or neg- 
ative, to the locomotion. For the three gaits considered here, the force is predominantly 
resistive, thus impeding locomotion. These computational results are in agreement with a 
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FIG. 6. (Color online) Distribution of the Tzz and r^^ components of the polymeric stress for pullers 
in the viscoelasic fluid with We = 0.5. Both the Tzz and Trz components lead to a net polymeric 
(elastic) force indicated by the gray arrows, hindering locomotion. The location of this elastic 
forces is right above the vortex ring, implying that they are generated by the elongational flow. 

recent experimental study of swimming Caenorhabditis elegans in polymeric fluids where the 
authors suggest that elongational viscosities can explain the measured decrease in swimming 
speeds [23] . 



V. SWIMMING POWER IN VISCOELASTIC FLUIDS 

The mechanical power, expended by the spherical squirmer is reported as a function 
of the Weissenberg number, We^ for the three different gaits under investigation in Fig. [7[ 
When considering dimensional quantities, we observe that both the puller and pusher require 
significantly more power than the neutral squirmer due to the strong vortex rings (see Fig.[T]). 
Indeed, as reported in Ref. [31j, V grows quadratically with the dipole magnitude a. In order 
to factor out the influence of the swimming gait on the power expenditure, thereby isolating 
the effect of elasticity on the energy budget, we normalize the power with its corresponding 
value in the Newtonian fluid, termed V^ew In Fig- we observe a significant reduction in 
power in the viscoelastic fluid for all gaits; V decreases first rapidly in the range We G [0, 2], 
and tends to an asymptotic value when further increasing the polymer relaxation time. The 
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FIG. 7. (Color online) Swimming power, for pusher (squares, green online), puller (triangles, 
blue online) and neutral (circles, red online) squirmers, nondimensionalized by the swimming power 
in the Newtonian case, T^New-The power consumption of swimming in the Newtonian fluid are equal 
to PNew = S/Stt (neutral squirmer) and PNew — 367r (both pusher and puller). 

limiting value of V for We 1 is close to half the value of the Newtonian case {We = 0), 
namely, 7^|vi^e>i/^New ^ 0.5. The result is in agreement with the asymptotic analysis of 
Ref. [2)J which predicts V\we:^i/T^New = 1^ W = 0.5 for the computations presented here). 

As we will show below, the component of the power consumption associated with the 
Newtonian solvent is almost constant with variations in We^ while the contribution of the 
polymeric stresses decreases significantly (up to one order of magnitude decrease at large 
We). Since, for /3 = 0.5, the Newtonian and polymeric contribution are the same for We = 0, 
the power at large Weissenberg numbers is about half that in a Newtonian fiuid. 

To probe the power saving for micros wimmers in viscoelastic fiuids, we decompose the 
power into three parts. First, we use the divergence theorem to transform the total power, 
7^, which is the integral of the stress at the surface times the swimming speed, into a volume 
integral as done in j7] 



where n is the unit normal outward the swimmer surface u is the velocity in the laboratory 




(9) 
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frame and cr is the stress tensor introduced above. Rearranging the integrand, one can write 

V • (u • r) = /5/xcc;2 + 2/3/x(Vu : Vu) + E : r^. (10) 

Consequently, V is found to be the sum of three contributions, V = Vq + Vdv + Vp\ a com- 
ponent related to the flow vorticity, Vq^ = Jy /3/xa;^rf]/ , a component related to the velocity 
gradient, Vdv = /y2;5/x(Vu : Vu)rfl/, and a polymeric component, Vp = /yE : r^dV^ 
where E is the rate-of-strain tensor. The component related to the velocity derivative, Vdv^ 
is only dependent on the tangential velocity distribution we impose on the surface and thus 
is only a function of a. This can be seen by re- writing 

Vdv = 2/3/x / (Vu) : (Vu) dV = -2/3/x / n • (u • Vu) dS. (11) 
Jv Js 




FIG. 8. (Color online) Contributions to the total power expended by swimming, versus Weis- 
senberg number We, /3 = 0.5, for neutral squirmer (left) and puller (right). The three contributions 
are defined as: Vq = Jy ^iioj'^dV, Vdv = /y 2/3/x(Vu : V\x)dV and Vp ^ E : r^dV. 

The three diflFerent contributions to the total consumed power are displayed in Fig. [8] as 
a function of the Weissenberg number. We, for two swimming gaits, neutral swimmer (left 
panel) and puller (right panel). Energetic results for pushers are qualitatively similar to 
those for pullers. First note that for We = and our choice /3 = 0.5, the power component 
due to the polymeric stress is Vp = Vq + Vdv- As seen from the deflnition of the polymeric 
model, reduces to Newtonian stress {I — f3 )(Vu + Vu^) when We = in Eq. ([s]). 

For a flxed mode of surface motion, Vdv^ is strictly constant with We since it depends 
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only on the values of the velocity at the surface (see Eq.ll). For both gaits, we observe 



computationally that the vorticity contribution, P^, shows little variations with We. The 
major contribution to the reduction in swimming power stems thus from the significant 
reduction of the polymeric part, Vp^ with increasing values of the Weissenberg number. 

The contribution of polymeric stresses to the consumed power is the double dot product 
of two tensors, the rate-of-strain tensor, E = |(Vu + (Vu)^), and the polymeric stress 
tensor, r^. Polymer chains get stretched by the rate-of-strain. In fiuids with low relaxation 
times, polymeric stretching occurs immediately, resulting in high spatial correlation between 
the polymeric stress tensor and the rate-of-strain tensor. Thus, the product of the two 
tensors takes large values. In contrast, for fiuids with large relaxation times, polymer chains 
get fully stretched at a position far away from where they get excited, and the spatial 
correlation of the two tensors decreases. Polymeric deformation and rate-of-strain become 
therefore separated in space as We increases, which gives a systematically smaller power 
consumption, as observed computationally. 



We = 0.1 



We 




FIG. 9. (Color online) Polymeric power density, Dp, for pushers in fluids with We = 0.1 (left) and 
We = 2 (right). Dotted elliptical circles mark areas with negative values oi Dp. The values of Dp 



along the swimmer surface are displayed in Fig. 10 



In Fig. [9| we plot the power density. Dp = 27r|r|E : r^, for the spherical pusher for 
two different values of We. When We = 0.1 (left half of the figure), large values oi Dp 
are observed in the region very close to the swimmer, especially in the equatorial region. 
This is the region where most of the energy is expended to elongate the polymers. A 
region of positive Dp also exists near the equatorial z = plane for the case with We = 2 
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shown on the right, although of reduced size. In addition, negative density can be seen 
on the planes z ^ ±0.4 (indicated by elliptical circles in the figure), which can also be 
clearly identified in Fig. [TO] where we display the value of Dp evaluated along the surface 
of the pusher. These areas of negative power density produce further decrease of the total 
expended power by providing energy to the system. This energy analysis allows us to 
understand the reduction in swimming power observed in viscoelastic fiuids for any type of 
surface deformation considered. 

We finally consider the situation of locomotion at constant power. Some bacteria have 
been observed to swim with kinematics consistent with constant-power conditions [50j , and 
it is thus interesting to consider the swimming speed which would be attained by squirmers 
if they were under the constraints of expending the same power as in a Newtonian fiuid (the 
results in Fig. [2] were obtained at constant intensity of the velocity at the boundary). The 
swimming speed obtained rescaling the previous data to obtain constant-power locomotion 



is displayed in Fig. [TT| To calculate the velocity at constant power, we need to increase the 
value of the surface velocity Bi to account for the power-saving in polymeric fiuid discussed 
above. If we use B[ to denote the new value of Si, constant power requirement is written 
as {B[) T^lv^gj^^/^ = BfV^ewi where We {B[) = = -^We] we use a simple interpolation 
from our computational results in order to obtain V\^^(^-Q,y Thus, the swimming speed with 



Dp 
2000 

1500 

1000 

500 



-500 



-We 


= 0.1! 

= 2 \ 




r/V) 






1 









TT 

4 



TT 

2 



4 



71 



e 



FIG. 10. (Color online) Spatial distribution of polymeric power density, Dp, along the surface of 
pushers in viscoelastic fluids with We = 0.1 (solid line, red online) and 2 (dashed line, blue online); 
6 is deflned in Sec. [IT} as the angle between the swimming direction and the position vector of each 
point. 
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FIG. 11. Swimming speed scaled to ensure locomotion at constant power, versus Weissenberg 
number, VFe, for /3 = 0.5. Nondimensional form is plotted for the three swimming gaits: 

neutral (circles, red online), puller (triangles, blue online) and pusher (squares, green online) 



U is found to 



constant power is computed as [/ = -^U\y^^^^,y As can be seen in Fig. 11 
be systematically larger than that in the Newtonian conditions for all the cases considered, 
and the swimming speed systematically increases with the Weissenberg number. The puller 
and the neutral squirmer are about 20% faster than the pusher. 



VI. SWIMMING EFFICIENCY IN VISCOELASTIC FLUIDS 



The hydrodynamic efficiency of squirmers in viscoelastic fluid is displayed in Fig. [T2j The 
efficiency is deflned as the ratio between the rate work needed to pull a sphere of same size 
in the same fluid at the swimming speed U and the swimming power 7^ of a self-propelled 
swimmer [5ll [52] 

FU , , 

V=^. (12) 

where F is the force required to drag the spherical squirmer body at the speed U. Similarly 
to what we did for the power consumption, the efficiency is normalized in Fig. [12] by its 
Newtonian value. We observe that all swimmers have a higher swimming efficiency than in 
a Newtonian fluid. For all three gaits, as We increases from 0, the scaled efficiency firstly 
increases rapidly with We^ then decreases asymptotically. We can thus identify an optimal 
We number corresponding to the maximum efficiency; this is around 2 for the neutral and 
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FIG. 12. (Color online) Swimming efficiency, 77, normalized by the corresponding value for New- 
tonian swimming, r/New? versus Weissenberg number. We. The efficiencies of swimming in the 
Newtonian fluid are equal to r/New — 0.5 (neutral squirmer) and r/New — 0.037 (both pusher and 
puller) . 

puller swimmer and around 0.5 for the pusher. 

VII. VELOCITY DECAY IN VISCOELASTIC FLUIDS 

Because of its implications to collective cell behavior [23l [5l] , it is of interest to investigate 
the effect of the fluid elasticity on the spatial decay of the flow perturbation induced by the 
swimming motion. We compute the axial velocity along the symmetry axis of the domain 
for the pusher and puller (similar analysis was carried out for the neutral swimmer in our 
previous paper [24j). The velocity decays with a power-law behavior, \u\ ^ l/r^, with 7 = 2 
in the Newtonian case for both pushers and pullers, as sufficiently far away from the cell the 
higher order terms (1/r^ and 1/r^) give negligible contribution to the flow field (7 = 3 in 
the neutral potential- flow case). In the non-Newtonian case, we estimate the value of 7 by 
fitting a power law to the numerical results in a measurement region extending from about 
r ^ 5D to the end of the computational domain. The estimated values of 7 are displayed in 
Fig. [13] as a function of the Weissenberg number. We. We see clearly that the velocity always 
decays faster in viscoelastic fluids than in the Newtonian fluid. For both pusher and puller, 
7 is not monotonically increasing when with We^ but instead reaches a maximum value at 
intermediate Weissenberg numbers. The maximum value of 7 is seen to take place as We is 
around 5 for the pusher and around 1 for the puller, which approximately corresponds with 
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FIG. 13. (Color online) Power law exponent, 7, for spatial decay of the axial velocity along the 
axis of the domain (r = 0) {\u\ ^ r~^, see text), as a function of the Weissenberg number. We. 
The inset plot highlights the difference between pusher and puller. 



the values of the Weissenberg numbers at which the swimming speed, [/, is minimum for 
both gaits (although slightly above). We thus observe a negative correlation between the 
swimming speed and the velocity decay rate. We further note that the velocity decay rate 
is larger for the pusher than for the puller, which is consistent with our observation of the 
negative correlation as the swimming speed of a pusher is smaller than that of a puller. 



VIII. STRESSLET IN VISCOELASTIC FLUIDS 

In this final section we address how the viscoelastic modification of the squirming mo- 
tion would affect the rheology of an active suspension of squirmers. When the size of the 
microorganisms is much smaller than the scale of the fiow field under consideration, it is 
convenient to use a continuum model of the active suspension. The first order correction to 
the bulk viscosity in terms of the volume fraction of swimmers is then given by the stresslet 
associated with an individual swimmer. Batchelor [53] derived a relation between the bulk 
stress and the conditions at the surfaces of individual particles embedded in a Newtonian 
solvent (his original formulation is given in dimensional form and we present the nondimen- 



sional version here using the scaling we introduced in Sec. II). The bulk stress a is given 
by 

<T = I.T. + 2E + (t', (13) 
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where I.T. is an isotropic tensor, i.e. I.T. = akk^^ E is the volume-averaged rate-of-strain 
tensor {Eij = ^ + ^.V^ where V is the volume occupied by the fluid and the 

particles) and represents the disturbance stress due to the presence of squirmers; see also 
Ref [54] for the application to active suspensions in Newtonian fluids. The particle bulk 
stress induced by force- and torque-free particles (such as swimmers) in the volume V 
can be expressed as 

-^ = ^ES, (14) 



V 

^ ^ dA. (15) 



- {(cr • n) X + X (cr • n)} — -X • cr • nl — (un + nu) 

'As [2 3 

The stress is thus given by the volume average of the so-called stresslet S, where the sum 



is taken over all the swimmers. The stresslet S is deflned Eq. (15) where a is the stress 
tensor, As is the surface of one of the squirmers, n is the unit outward normal vector, x is 
the position vector and I is the unit tensor. 

Ishikawa and coworkers [17j exploited the formulation above to study the stresslet S 
of a single squirmer and a pair of squirmers in the Newtonian fluid both analytically and 
numerically. They showed that for single squirmers 

S = ^7r(3ee-I)52, (16) 

where e is the swimming direction and B2 the second squirming mode deflned in SecjlT} The 
assumption of axisymmetric bodies makes the off-diagonal part of S zero. In a Newtonian 
fluid, the flrst mode Bi determines the swimming speed while B2 the stresslet magnitude; 
the stresslet of a neutral squirmer {B2 = 0) is therefore zero. 

In order to derive a similar formalism for a squirmer in a non-Newtonian flow, we 
consider the additional polymeric contribution, r^, to the Newtonian stress, a = —pi + 
/3 (Vu + Vu^) + tP. The same derivation as in Refs. |53l [55] can be performed for non- 
Newtonian fluids, see for example Ref. [56j. In our case, we seek an expression for 

6- = I.T. + 2/3E + rp + o■^ (17) 

where fp is the volume average of the polymeric stress (r^ = ^ TpdV) . The total bulk 



20 



stress can be written as 



^ [ {-pI + /5(Vu + Vu^)+rJdy + ^ / adV, 

y Jv-Vs ^ Jvs 



(18) 



where Vs is the volume occupied by the squirmer. Combining Eq. (17) and Eq. (18), we 
have 



^ JVq ^ JVq ^ JV 



(19) 

'Vs JVs JVs 

The volume integral of the Newtonian contribution to the total stress inside the squirmer 
is re-written as a surface integral by applying the divergence theorem: the first term on the 



right hand side of Eq. (19) becomes J^^ (5 (un + nu) dA^ whereas the total stress inside 
the particles contributes to the stresslet 



f adV= [ (o- n)xrfA- / (V • o-) xrfy, 

Jvs JAs JVs 



(20) 



where we have assumed V • cr = inside the solid body. For a swimming squirmer which is 
force-free and torque-free, (cr • n) x = | {(cr • n) x + x (cr • n)} — |x • cr • nl as in Ref. [TTj. 

Finally we obtain the new definition of the stresslet 



As 



^ {(o- • n) X + X (o- • n)} - ^x • o- • nl 



(3 (un + nu) 



dA 



T^dV. (21) 



Vs 



The stresslets S computed from Eq. ( [21] ) for the pusher {a = —5), puller {a = 5) and 

Because of axisymmetry, there are only 
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neutral squirmer {a = 0) are displayed in Fig. 
two non-zero components of S, namely, the rr component in the direction normal to the 
swimming direction and the zz component in the swimming direction. 

First, we see that for the pusher and puller the magnitude of both the rr and zz compo- 
nents decreases monotonically with increasing We^ and reaches an asymptotic value (while 
keeping the same sign). This corresponds to a maximum decrease of approximately 20% 
with respect to the Newtonian value at We = 10. This implies that a dilute suspension of 
swimming cells will have a weaker rheological effect in a viscoelastic fiuid. 

The decrease in the stresslet can tentatively be explained by the fact that the polymeric 
stress always opposes the fiow which generates it. Let us consider the zz component of a 
pusher as an example; in that case, the tangential velocity on the surface pushes the fiuid 
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away from the cell along the swimming direction as shown in Fig. [T} resulting in an exten- 
sional flow in front and behind the cell body with a consequent negative zz component of the 
stresslet. Simultaneously, the extensional flow produces strong positive polymeric stresses 
(see Fig. [i]) counteracting the eflFect created by the wall actuation, and hence decreasing the 
net stress; the impact on the rr component can be explained in a similar fashion. Note that 
for both the pusher and puller, the magnitude of the zz component is approximatively twice 



that of the rr component, as can be seen from Eq. (16). 



Interestingly, we also observe in Fig. [14] that for the neutral squirmer, which has zero 
stresslet in the Newtonian limit, viscoelastic stresses induce a non-zero stresslet. The value of 
the rr component is lower than that of the zz component and of same sign. Considering two 
squirmers swimming in the same direction side by side, at leading order in their separation 
distance they have no influence on each other in a Newtonian fluid, but will repel each 
other due to the negative value of the rr component of the stresslet. Similarly, two neutral 
squirmers swimming in line along the same swimming direction will tend to increase their 
relative distance owing to the negative value of the zz component of the stresslet. The 
peak value of the rr and zz components occur at We ^ 1, at which the minimum swimming 
speed of the neutral squirmer appears, as well as the fastest velocity decay. Interestingly, this 
viscoelastic repulsion between swimming squirmers has the sign opposite to the viscoelastic 
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FIG. 14. (Color online) Stresslets induced by swimmers in viscoelastic fluids: rr and zz components 



of the tensor S, Eq. (21), for different swimming gaits. Left: pusher (a — —5) and puller (a — 5). 
Right: neutral squirmer (o: = 0). 
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attraction observed between sedimenting spheres in polymeric fluids. 



IX. CONCLUSION 

In the current paper, we numerically compute the swimming speed, power, and hydrody- 
namic efliciencies of model cihated cells swimming in a viscoelastic Giesekus fluid. We use 
the squirmer model where the self-propulsion is driven by tangential surface deformation 
and consider three types of boundary conditions to mimic swimming gaits typical of cell 
locomotion (potential squirmer, pusher, and puller) all having the same swimming speed in 
a Newtonian fluid. The characteristics of the flow and of the polymeric stress are examined 
to help understand our computational results. 

For constant magnitude of the surface deformation velocity, we show that the swimming 
speed decreases in the viscoelastic fluid compared its the Newtonian value for all cases 
considered, with a minimum obtained at intermediate values of the Weissenberg number. 
Neutral swimmer and puller have the highest velocity, whereas the pusher is about 15% 
slower. For the potential squirmer and the pusher, the reduced self-propulsion is explained 
by the accumulation of highly stretched polymers behind the body. For low values of the 
polymer relaxation times. We < 1, the stretching of the polymers increases hence the initial 
decrease in the swimming speed. For the largest values of the Weissenberg number, however, 
the region of stretched polymers behind the swimmer becomes thinner and the induced 
elastic resistance decreases, thus explaining the slow recovery of the swimming speed. In 
the case of pullers, a different velocity field is created and stretched polymers are advected 
away from the cell on the sides. This induces a force with a component both in the radial 
and axial direction, the latter explaining the reduced swimming speed. Similar kinematic 
analysis could be used to provide detailed guidelines for the design of efficient swimmers 
in viscoelastic ffuids - the basic design principle consisting in reducing the accumulation of 
stretched polymers at the cell surface in the direction of motion. 

Along with the decrease in swimming speed, we observe a significant reduction of swim- 
ming power in all cases. To understand this observation, we analytically decompose the 
expression for the power expended by cells swimming in viscoelastic fluids into three con- 
tributions (integrals over the flow domain): one associated to the vorticity induced in the 
flow, one only set by the surface deformation, and one related to the polymeric stresses. In 
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this way, we are able to clearly identify the reduction in consumed power with a reduction 
in the polymeric contribution and interpret the observed reduction in power as due to an 
increasing spatial de-correlation between the flow shear and the induced polymer deforma- 
tion. In addition, we observe regions of negative power density indicating that the polymers 
in the fluid have the capacity to flrst transform the mechanical energy into potential energy 
(polymer deformation), and then feed it back to the fluid thus reducing the consumed power. 
We also consider the hydrodynamic efficiencies and show that self-propulsion is always more 
efficient in a viscoelastic fluid, with a maximum relative efficiency obtained for intermediate 
values of the Weissenberg number. 

We further investigate the influence of the fluid elasticity on the decay rate of the pertur- 
bation velocity induced by the squirmer. We flnd that the decay rate is faster in viscoelastic 
fluid and is larger for the pusher than for the puller. Finally, to address the influence of 
viscoelasticity on the rheology of an active suspension of swimmers, we extract from the 
computational results the stresslet associated with an isolated swimmer. Notably, viscoelas- 
ticity induces a non-zero stresslet for the potential squirmer (rigorously equal to zero in 
the Newtonian case). As a difference, the magnitude of the stresslet for both pushers and 
pullers decrease by about 20%. These results presents a flrst quantitative step toward an 
understanding of the dynamics and rheology of active suspensions in viscoelastic fluids. 
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